rm(list = ls())

install.packages(c("gsynth", "remotes", "did", "haven", "panelView" , "Rcpp", "ggplot2", "GGally", "foreach", "doParallel", "abind", "PanelMatch", "dotwhisker", "broom", "dplyr", "metafor", "foreign", "estimatr", "purrr", "brms", "tidybayes", "tidyr"))
install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
cmdstanr::install_cmdstan()

library(cmdstanr)
library(remotes)
library(haven)
library(did)
library(gsynth)
library(panelView)
require(Rcpp) 
require(ggplot2)  
require(GGally) 
require(foreach)  
require(doParallel) 
require(abind) 
library(PanelMatch)
library(dotwhisker)
library(broom)
library(dplyr)
library(metafor)
library(foreign)
library(estimatr)
library(purrr)
library(brms)
library(tidybayes)
library(tidyr)
set.seed(8675309)

setwd("/Users/cb2257/Desktop/Replication/")

Study1 <- read.csv(file="data/Study1Data.csv")
Study2 <- read.csv(file="data/Study2Data.csv")
Study3 <- read.csv(file="data/Study3Data.csv")
options(scipen=999)

###Estimate Effect of Analogies Relative to the Control (Study 1)

Study1$MunichAnalogy <- NA 
Study1$MunichAnalogy[Study1$Treatment=="Control (Munich)"] <- 0 
Study1$MunichAnalogy[Study1$Treatment=="Treatment (Munich)"] <- 1 

Study1$CubanMissileAnalogy <- NA 
Study1$CubanMissileAnalogy[Study1$Treatment=="Control (Cuban)"] <- 0 
Study1$CubanMissileAnalogy[Study1$Treatment=="Treatment (Cuban)"] <- 1 

Study1$success_pooled_munich <- (Study1$General_S + Study1$Tactical_S_M + Study1$Strategic) / 3

Study1$success_pooled_cuban <- (Study1$General_S + Study1$Tactical_S_C + Study1$Strategic) / 3

Study1$success_pooled_munich_binary <- (Study1$General_S_Binary + Study1$Tactical_S_M_Binary + Study1$Strategic_Binary) /3

Study1$success_pooled_cuban_binary <- (Study1$General_S_Binary + Study1$Tactical_S_C_Binary + Study1$Strategic_Binary) /3

Study1$CB_Pooled_Binary <- (Study1$CB_CostsLow_Binary + Study1$CB_BenefitsHigh_Binary + Study1$CB_BenefitsOutweighCosts_Binary) / 3

Study1$Conf_Pooled_Binary <- (Study1$Conf_KnowHistory_Binary + Study1$Conf_Intelligent_Binary + Study1$Conf_Competent_Binary) / 3

model1_cred <- lm_robust(Cred ~ MunichAnalogy, data=Study1)

model2_cred <- lm_robust(Cred ~ CubanMissileAnalogy, data=Study1)

model1_cred_binary <- lm_robust(Cred_Binary ~ MunichAnalogy, data=Study1)

model2_cred_binary <- lm_robust(Cred_Binary ~ CubanMissileAnalogy, data=Study1)

model1_success <- lm_robust(success_pooled_munich ~ MunichAnalogy, data=Study1)

model2_success <- lm_robust(success_pooled_cuban ~ CubanMissileAnalogy, data=Study1)

model1_success_binary <- lm_robust(success_pooled_munich_binary ~ MunichAnalogy, data=Study1)

model2_success_binary <- lm_robust(success_pooled_cuban_binary ~ CubanMissileAnalogy, data=Study1)

model1_cb <- lm_robust(CB_Pooled ~ MunichAnalogy, data=Study1)

model2_cb <- lm_robust(CB_Pooled ~ CubanMissileAnalogy, data=Study1)

model1_cb_binary <- lm_robust(CB_Pooled_Binary ~ MunichAnalogy, data=Study1)

model2_cb_binary <- lm_robust(CB_Pooled_Binary ~ CubanMissileAnalogy, data=Study1)

model1_traits <- lm_robust(Conf_Pooled ~ MunichAnalogy, data=Study1)

model2_traits <- lm_robust(Conf_Pooled ~ CubanMissileAnalogy, data=Study1)

model1_traits_binary <- lm_robust(Conf_Pooled_Binary ~ MunichAnalogy, data=Study1)

model2_traits_binary <- lm_robust(Conf_Pooled_Binary ~ CubanMissileAnalogy, data=Study1)

model1_moral <- lm_robust(M_M ~ MunichAnalogy, data=Study1)

model2_moral <- lm_robust(M_C ~ CubanMissileAnalogy, data=Study1)

model1_moral_binary <- lm_robust(M_M_Binary ~ MunichAnalogy, data=Study1)

model2_moral_binary <- lm_robust(M_C_Binary ~ CubanMissileAnalogy, data=Study1)

###Estimate Effect of Analogies Relative to the Control (Study 2)

Study2$BosniaAnalogy <- NA 
Study2$BosniaAnalogy[Study2$Z1==0] <- 0 
Study2$BosniaAnalogy[Study2$Z1==1 | Study2$Z1==2] <- 1 

Study2$VietnamAnalogy <- NA 
Study2$VietnamAnalogy[Study2$Z1==3] <- 0 
Study2$VietnamAnalogy[Study2$Z1==4] <- 1 

Study2$PhilippinesAnalogy <- NA 
Study2$PhilippinesAnalogy[Study2$Z1==3] <- 0 
Study2$PhilippinesAnalogy[Study2$Z1==5] <- 1 

Study2$General_S <- as.numeric(Study2$General_S)

Study2$Tactical_Success_Valence <- as.numeric(Study2$Tactical_Success_Valence)

Study2$Strategic_Success_Valence <- as.numeric(Study2$Strategic_Success_Valence)

Study2$success_pooled_bosnia <- (Study2$General_S + Study2$Tactical_Success_Valence + Study2$Strategic_Success_Valence) / 3

Study2$Tactical_Familiar <- as.numeric(Study2$Tactical_Familiar)

Study2$Strategic_Success_Familiar <- as.numeric(Study2$Strategic_Success_Familiar)

Study2$success_pooled_iran_scenario <- (Study2$General_S + Study2$Tactical_Familiar + Study2$Strategic_Success_Familiar) / 3

Study2$General_S_Binary <- as.numeric(Study2$General_S_Bi.ry)

Study2$Tactical_Success_Valence_Binary <- as.numeric(Study2$Tactical_Success_Valence_Bi.ry)

Study2$Strategic_Success_Valence_Binary <- as.numeric(Study2$Strategic_Success_Valence_Bi.ry)

Study2$success_pooled_bosnia_binary <- (Study2$General_S_Binary + Study2$Tactical_Success_Valence_Binary + Study2$Strategic_Success_Valence_Binary) / 3

Study2$Tactical_Familiar_Binary <- as.numeric(Study2$Tactical_Familiar_Bi.ry)

Study2$Strategic_Success_Familiar_Binary <- as.numeric(Study2$Strategic_Success_Familiar_Bi.ry)

Study2$success_pooled_iran_scenario_binary <- (Study2$General_S_Binary + Study2$Tactical_Familiar_Binary + Study2$Strategic_Success_Familiar_Binary) / 3

Study2$CB_BenefitsHigh_Bi.ry <- as.numeric(Study2$CB_BenefitsHigh_Bi.ry)

Study2$CB_BenefitsOutweighCosts_Bi.ry <- as.numeric(Study2$CB_BenefitsOutweighCosts_Bi.ry )

Study2$CB_CostsLow_Bi.ry <- as.numeric(Study2$CB_CostsLow_Bi.ry)

Study2$CB_Pooled_Binary <- (Study2$CB_BenefitsHigh_Bi.ry + Study2$CB_BenefitsOutweighCosts_Bi.ry + Study2$CB_CostsLow_Bi.ry) /3 

Study2$Conf_Competent_Bi.ry <- as.numeric(Study2$Conf_Competent_Bi.ry)
Study2$Conf_Intelligent_Bi.ry <- as.numeric(Study2$Conf_Intelligent_Bi.ry)
Study2$Conf_KnowHistory_Bi.ry <- as.numeric(Study2$Conf_KnowHistory_Bi.ry)

Study2$Conf_Pooled_Binary <- (Study2$Conf_Competent_Bi.ry + Study2$Conf_Intelligent_Bi.ry + Study2$Conf_KnowHistory_Bi.ry) / 3

model3_cred <- lm_robust(Cred_Val ~ BosniaAnalogy, data=Study2)

model3_cred_binary <- lm_robust(Cred__Val_Bi.ry ~ BosniaAnalogy, data=Study2)

Study2$Cred_Fam <- as.numeric(Study2$Cred_Fam)

model4_cred <- lm_robust(Cred_Fam ~ VietnamAnalogy, data=Study2)

model5_cred <- lm_robust(Cred_Fam ~ PhilippinesAnalogy, data=Study2)

Study2$Cred__Fam_Bi.ry <- as.numeric(Study2$Cred__Fam_Bi.ry)

model4_cred_binary <- lm_robust(Cred__Fam_Bi.ry ~ VietnamAnalogy, data=Study2)

model5_cred_binary <- lm_robust(Cred__Fam_Bi.ry ~ PhilippinesAnalogy, data=Study2)

model3_success <- lm_robust(success_pooled_bosnia ~ BosniaAnalogy, data=Study2)

model4_success <- lm_robust(success_pooled_iran_scenario ~ VietnamAnalogy, data=Study2)

model5_success <- lm_robust(success_pooled_iran_scenario ~ PhilippinesAnalogy, data=Study2)

model3_success_binary <- lm_robust(success_pooled_bosnia_binary ~ BosniaAnalogy, data=Study2)

model4_success_binary <- lm_robust(success_pooled_iran_scenario_binary ~ VietnamAnalogy, data=Study2)

model5_success_binary <- lm_robust(success_pooled_iran_scenario_binary ~ PhilippinesAnalogy, data=Study2)

Study2$CB_Pooled <- as.numeric(Study2$CB_Pooled)

model3_cb <- lm_robust(CB_Pooled ~ BosniaAnalogy, data=Study2)

model4_cb <- lm_robust(CB_Pooled ~ VietnamAnalogy, data=Study2)

model5_cb <- lm_robust(CB_Pooled ~ PhilippinesAnalogy, data=Study2)

model3_cb_binary <- lm_robust(CB_Pooled_Binary ~ BosniaAnalogy, data=Study2)

model4_cb_binary <- lm_robust(CB_Pooled_Binary ~ VietnamAnalogy, data=Study2)

model5_cb_binary <- lm_robust(CB_Pooled_Binary ~ PhilippinesAnalogy, data=Study2)

Study2$Conf_Pooled <- as.numeric(Study2$Conf_Pooled)

model3_traits <- lm_robust(Conf_Pooled ~ BosniaAnalogy, data=Study2)

model4_traits <- lm_robust(Conf_Pooled ~ VietnamAnalogy, data=Study2)

model5_traits <- lm_robust(Conf_Pooled ~ PhilippinesAnalogy, data=Study2)

model3_traits_binary <- lm_robust(Conf_Pooled_Binary ~ BosniaAnalogy, data=Study2)

model4_traits_binary <- lm_robust(Conf_Pooled_Binary ~ VietnamAnalogy, data=Study2)

model5_traits_binary <- lm_robust(Conf_Pooled_Binary ~ PhilippinesAnalogy, data=Study2)

Study2$Moral_Valence <- as.numeric (Study2$Moral_Valence)

Study2$Moral_Valence_Bi.ry <- as.numeric(Study2$Moral_Valence_Bi.ry)

model3_moral <- lm_robust(Moral_Valence ~ BosniaAnalogy, data=Study2)

model3_moral_binary <- lm_robust(Moral_Valence_Bi.ry ~ BosniaAnalogy, data=Study2)

Study2$Moral_Familiar <- as.numeric (Study2$Moral_Familiar)

Study2$Moral_Familiar_Bi.ry <- as.numeric(Study2$Moral_Familiar_Bi.ry)

model4_moral <- lm_robust(Moral_Familiar ~ VietnamAnalogy, data=Study2)

model5_moral <- lm_robust(Moral_Familiar ~ PhilippinesAnalogy, data=Study2)

model4_moral_binary <- lm_robust(Moral_Familiar_Bi.ry ~ VietnamAnalogy, data=Study2)

model5_moral_binary <- lm_robust(Moral_Familiar_Bi.ry ~ PhilippinesAnalogy, data=Study2)


###Estimate Effect of Analogies Relative to the Control (Study 3)

Study3$MunichAnalogy2 <- NA 
Study3$MunichAnalogy2[Study3$Treatment=="Control Taiwan"] <- 0 
Study3$MunichAnalogy2[Study3$Treatment=="Treatment Munich"] <- 1 

Study3$Cred <- as.numeric(Study3$Cred)

model6_cred <- lm_robust(Cred ~ MunichAnalogy2, data=Study3)

Study3$Cred_Binary <- NA
Study3$Cred_Binary[Study3$Cred==1 | Study3$Cred==2 | Study3$Cred==3] <- 0
Study3$Cred_Binary[Study3$Cred==4 | Study3$Cred==5] <- 1

model6_cred_binary <- lm_robust(Cred_Binary ~ MunichAnalogy2, data=Study3)

Study3$General_S <- as.numeric(Study3$General_S)

Study3$Tactical_S_M <- as.numeric(Study3$Tactical_S_M)

Study3$Strategic <- as.numeric(Study3$Strategic)

Study3$success_pooled <- (Study3$General_S + Study3$Tactical_S_M + Study3$Strategic) / 3

Study3$General_S_Bi.ry <- as.numeric(Study3$General_S_Bi.ry)

Study3$Tactical_S_M_Bi.ry <- as.numeric(Study3$Tactical_S_M_Bi.ry)

Study3$Strategic_Bi.ry <- as.numeric(Study3$Strategic_Bi.ry)

Study3$success_pooled_binary <- (Study3$General_S_Bi.ry + Study3$Tactical_S_M_Bi.ry + Study3$Strategic_Bi.ry) / 3

model6_success <- lm_robust(success_pooled ~ MunichAnalogy2, data=Study3)

model6_success_binary <- lm_robust(success_pooled_binary ~ MunichAnalogy2, data=Study3)

Study3$CB_Pooled <- as.numeric(Study3$CB_Pooled)

Study3$CB_BenefitsHigh_Bi.ry <- as.numeric(Study3$CB_BenefitsHigh_Bi.ry)

Study3$CB_BenefitsOutweighCosts_Bi.ry <- as.numeric(Study3$CB_BenefitsOutweighCosts_Bi.ry)

Study3$CB_CostsLow_Bi.ry <- as.numeric(Study3$CB_CostsLow_Bi.ry)

Study3$CB_Pooled_Binary <- (Study3$CB_BenefitsHigh_Bi.ry + Study3$CB_BenefitsOutweighCosts_Bi.ry + Study3$CB_CostsLow_Bi.ry) / 3

model6_cb <- lm_robust(CB_Pooled ~ MunichAnalogy2, data=Study3)

model6_cb_binary <- lm_robust(CB_Pooled_Binary ~ MunichAnalogy2, data=Study3)

Study3$Conf_Pooled <- as.numeric(Study3$Conf_Pooled)

Study3$Conf_Competent_Bi.ry <- as.numeric(Study3$Conf_Competent_Bi.ry)

Study3$Conf_Intelligent_Bi.ry <- as.numeric(Study3$Conf_Intelligent_Bi.ry)

Study3$Conf_KnowHistory_Bi.ry <- as.numeric(Study3$Conf_KnowHistory_Bi.ry)

Study3$Conf_Pooled_Binary <- (Study3$Conf_Competent_Bi.ry + Study3$Conf_Intelligent_Bi.ry + Study3$Conf_KnowHistory_Bi.ry) / 3

model6_traits <- lm_robust(Conf_Pooled ~ MunichAnalogy2, data=Study3)

model6_traits_binary <- lm_robust(Conf_Pooled_Binary ~ MunichAnalogy2, data=Study3)

Study3$M_M <- as.numeric(Study3$M_M)

Study3$M_M_Binary <- NA 
Study3$M_M_Binary[Study3$M_M==1 | Study3$M_M==2 | Study3$M_M==3] <- 0
Study3$M_M_Binary[Study3$M_M==4 | Study3$M_M==5] <- 1

model6_moral <- lm_robust(M_M ~ MunichAnalogy2, data=Study3)

model6_moral_binary <- lm_robust(M_M_Binary ~ MunichAnalogy2, data=Study3)

###Meta-Analysis Credibility 

meta_df <- 
  list(model1_cred, model2_cred, model3_cred, model4_cred, model5_cred, model6_cred) %>% 
  map_df(tidy) %>% 
  filter(term == "MunichAnalogy" | term == "CubanMissileAnalogy" | term =="BosniaAnalogy" |
           term == "VietnamAnalogy" | term == "PhilippinesAnalogy" | term == "MunichAnalogy2")

#Meta Effect = 0.197 and p-value = 0.0005
meta_fit <- rma.uni(yi = estimate, sei = std.error, method = "REML", data = meta_df)

ggdf <- meta_df %>% 
  bind_rows(tidy(meta_fit, conf.int = TRUE)) %>% 
  mutate(model = factor(c("Munich Analogy (Study 1)",
                          "Cuba Missile Crisis Analogy (Study 1)",
                          "Bosnian War Analogy (Study 2)",
                          "Vietnam War Analogy (Study 2)",
                          "Philippines War Analogy (Study 2)",
                          "Munich Analogy (Study 3)",
                          "Meta Average"),
                        levels = c("Munich Analogy (Study 3)",
                                   "Philippines War Analogy (Study 2)",
                                   "Vietnam War Analogy (Study 2)",
                                   "Bosnian War Analogy (Study 2)",
                                   "Cuba Missile Crisis Analogy (Study 1)",
                                   "Munich Analogy (Study 1)",
                                   "Meta Average")),
         Number = rep(1, 7))

# Assuming ggdf is your data frame
ggdf$is_meta <- ifelse(ggdf$model == "Meta Average", TRUE, FALSE)

plot <- ggplot(ggdf, aes(estimate, model)) +
  geom_point(aes(color = is_meta, shape = is_meta), size = 4) +
  geom_linerange(aes(xmin = conf.low, xmax = conf.high, color = is_meta), size = 1) +
  scale_color_manual(values = c("TRUE" = "black", "FALSE" = "black")) +
  scale_shape_manual(values = c("TRUE" = 17, "FALSE" = 16)) +  # Different shape for meta
  theme_minimal(base_size = 15) +
  ylab("") +
  xlab("\n Effect of Historical Analogies on \n Perceptions the President Chose the \n Best Foreign Policy Strategy") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red", size = 1) +
  geom_hline(yintercept = 6.5, linetype = "dotted", color = "black", size = 1) +
  theme(
    axis.text.x = element_text(size = 15, color = "black"),
    axis.text.y = element_text(size = 15, color = "black"),
    axis.title.y = element_text(size = 18, face = "bold"),
    axis.title.x = element_text(size = 18, face = "bold"),
    plot.caption = element_text(size = 15),
    legend.position = "none",  # Remove legend if not needed
    panel.background = element_rect(fill = "white", color = "white"),
    plot.background = element_rect(fill = "white", color = "white")
  ) +
  labs(
    title = ""
  ) +
  theme(axis.text.y = element_text(face = ifelse(levels(ggdf$model) == "Meta Average", "bold", "plain")))

# Display the plot
print(plot)

ggsave("figures/effect_of_historical_analogies.png", plot = plot, width = 10, height = 7, dpi = 300, device = "png")


###Meta-Analysis Credibility BINARY

meta_df <- 
  list(model1_cred_binary, model2_cred_binary, model3_cred_binary, model4_cred_binary, model5_cred_binary, model6_cred_binary) %>% 
  map_df(tidy) %>% 
  filter(term == "MunichAnalogy" | term == "CubanMissileAnalogy" | term =="BosniaAnalogy" |
           term == "VietnamAnalogy" | term == "PhilippinesAnalogy" | term == "MunichAnalogy2")

#Meta Effect = 0.074 and p-value = 0.0005
meta_fit <- rma.uni(yi = estimate, sei = std.error, method = "REML", data = meta_df)

ggdf <- meta_df %>% 
  bind_rows(tidy(meta_fit, conf.int = TRUE)) %>% 
  mutate(model = factor(c("Munich Analogy (Study 1)",
                          "Cuba Missile Crisis Analogy (Study 1)",
                          "Bosnian War Analogy (Study 2)",
                          "Vietnam War Analogy (Study 2)",
                          "Philippines War Analogy (Study 2)",
                          "Munich Analogy (Study 3)",
                          "Meta"),
                        levels = c("Munich Analogy (Study 3)",
                                   "Philippines War Analogy (Study 2)",
                                   "Vietnam War Analogy (Study 2)",
                                   "Bosnian War Analogy (Study 2)",
                                   "Cuba Missile Crisis Analogy (Study 1)",
                                   "Munich Analogy (Study 1)",
                                   "Meta")),
         Number = rep(1,7))


# Multiply x-axis values (estimate) by 100
ggdf$is_meta <- ifelse(ggdf$model == "Meta Average", TRUE, FALSE)

ggdf$estimate <- ggdf$estimate * 100

# Plot with modified x-axis
plot <- ggplot(ggdf, aes(estimate, model)) +
  geom_point(aes(color = is_meta, shape = is_meta), size = 4) +
  geom_linerange(aes(xmin = conf.low * 100, xmax = conf.high * 100, color = is_meta), size = 1) +
  scale_color_manual(values = c("TRUE" = "black", "FALSE" = "black")) +
  scale_shape_manual(values = c("TRUE" = 17, "FALSE" = 16)) +  # Different shape for meta
  theme_minimal(base_size = 15) +
  ylab("") +
  xlab("\n Effect of Historical Analogies on \n Perceptions the President Chose the \n Best Foreign Policy Strategy") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red", size = 1) +
  geom_hline(yintercept = 6.5, linetype = "dotted", color = "black", size = 1) +
  theme(
    axis.text.x = element_text(size = 15, color = "black"),
    axis.text.y = element_text(size = 15, color = "black"),
    axis.title.y = element_text(size = 18, face = "bold"),
    axis.title.x = element_text(size = 18, face = "bold"),
    plot.caption = element_text(size = 15),
    legend.position = "none"  # Remove legend if not needed
  ) +
  labs(
    title = ""
  ) +
  theme(axis.text.y = element_text(face = ifelse(levels(ggdf$model) == "Meta Average", "bold", "plain")))

# Display the plot
print(plot)


###Meta-Analysis Mechanisms 

##Success
meta_df <- 
  list(model1_success, model2_success, model3_success, model4_success, model5_success, model6_success) %>% 
  map_df(tidy) %>% 
  filter(term == "MunichAnalogy" | term == "CubanMissileAnalogy" | term =="BosniaAnalogy" |
           term == "VietnamAnalogy" | term == "PhilippinesAnalogy" | term == "MunichAnalogy2")

#Meta Effect = 0.192 and p-value = <0.0001
meta_success <- rma.uni(yi = estimate, sei = std.error, method = "REML", data = meta_df)

##Cost-Benefit
meta_df <- 
  list(model1_cb, model2_cb, model3_cb, model4_cb, model5_cb, model6_cb) %>% 
  map_df(tidy) %>% 
  filter(term == "MunichAnalogy" | term == "CubanMissileAnalogy" | term =="BosniaAnalogy" |
           term == "VietnamAnalogy" | term == "PhilippinesAnalogy" | term == "MunichAnalogy2")

#Meta Effect = 0.138 and p-value = <0.0001
meta_cb <- rma.uni(yi = estimate, sei = std.error, method = "REML", data = meta_df)


##Presidential Traits
meta_df <- 
  list(model1_traits, model2_traits, model3_traits, model4_traits, model5_traits, model6_traits) %>% 
  map_df(tidy) %>% 
  filter(term == "MunichAnalogy" | term == "CubanMissileAnalogy" | term =="BosniaAnalogy" |
           term == "VietnamAnalogy" | term == "PhilippinesAnalogy" | term == "MunichAnalogy2")

#Meta Effect = 0.299 and p-value = <0.0001
meta_traits <- rma.uni(yi = estimate, sei = std.error, method = "REML", data = meta_df)


##Moral Obligation
meta_df <- 
  list(model1_moral, model2_moral, model3_moral, model4_moral, model5_moral, model6_moral) %>% 
  map_df(tidy) %>% 
  filter(term == "MunichAnalogy" | term == "CubanMissileAnalogy" | term =="BosniaAnalogy" |
           term == "VietnamAnalogy" | term == "PhilippinesAnalogy" | term == "MunichAnalogy2")

#Meta Effect = 0.117 and p-value = <0.0001
meta_moral <- rma.uni(yi = estimate, sei = std.error, method = "REML", data = meta_df)

meta_results <- data.frame(
  estimate = c(meta_success$beta, meta_cb$beta, meta_traits$beta, meta_moral$beta),
  conf.low = c(meta_success$ci.lb, meta_cb$ci.lb, meta_traits$ci.lb, meta_moral$ci.lb),
  conf.high = c(meta_success$ci.ub, meta_cb$ci.ub, meta_traits$ci.ub, meta_moral$ci.ub),
  category = c("Likelihood of Policy Success", "Policy's Benefits Outweigh its Costs", "Positive Presidential Traits", "Moral Obligation to Intervene")
)

# Order the categories for plotting
meta_results$category <- factor(meta_results$category, levels = c("Moral Obligation to Intervene", "Positive Presidential Traits", "Policy's Benefits Outweigh its Costs", "Likelihood of Policy Success"))

plot <- ggplot(meta_results, aes(x = estimate, y = category, color = category)) +
  geom_point(size = 4) +
  geom_linerange(aes(xmin = conf.low, xmax = conf.high), size = 1) +
  theme_minimal(base_size = 15) +
  ylab("") +
  xlab("\n Effect of Historical Analogies") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red", size = 1) +
  theme(
    axis.text.x = element_text(size = 15, color = "black"),
    axis.text.y = element_text(size = 15, color = "black"),
    axis.title.y = element_text(size = 18, face = "bold"),
    axis.title.x = element_text(size = 18, face = "bold"),
    plot.caption = element_text(size = 15),
    panel.background = element_rect(fill = "white", color = "white"),
    plot.background = element_rect(fill = "white", color = "white")
  ) +
  labs(
    title = "",
  ) +
  scale_color_manual(values = c("Likelihood of Policy Success" = "black", "Policy's Benefits Outweigh its Costs" = "black", "Positive Presidential Traits" = "black", "Moral Obligation to Intervene" = "black")) +
  theme(legend.position = "none")+ scale_x_continuous(
  breaks = c(-0.2, -0.2, 0, 0.2, 0.4),  # whatever tick marks you want
  limits = c(-0.1, 0.5)                 # full axis range
)

# Display the plot
print(plot)

ggsave("figures/mechanisms.png", plot = plot, width = 10, height = 7, dpi = 300, device = "png")

###Meta-Analysis Mechanisms (BINARY)

##Success
meta_df <- 
  list(model1_success_binary, model2_success_binary, model3_success_binary, model4_success_binary, model5_success_binary, model6_success_binary) %>% 
  map_df(tidy) %>% 
  filter(term == "MunichAnalogy" | term == "CubanMissileAnalogy" | term =="BosniaAnalogy" |
           term == "VietnamAnalogy" | term == "PhilippinesAnalogy" | term == "MunichAnalogy2")

#Meta Effect = 0.0775 and p-value = <0.0001
meta_success <- rma.uni(yi = estimate, sei = std.error, method = "REML", data = meta_df)


##Cost-Benefit
meta_df <- 
  list(model1_cb_binary, model2_cb_binary, model3_cb_binary, model4_cb_binary, model5_cb_binary, model6_cb_binary) %>% 
  map_df(tidy) %>% 
  filter(term == "MunichAnalogy" | term == "CubanMissileAnalogy" | term =="BosniaAnalogy" |
           term == "VietnamAnalogy" | term == "PhilippinesAnalogy" | term == "MunichAnalogy2")

#Meta Effect = 0.0581 and p-value = <0.0001
meta_cb <- rma.uni(yi = estimate, sei = std.error, method = "REML", data = meta_df)


##Presidential Traits
meta_df <- 
  list(model1_traits_binary, model2_traits_binary, model3_traits_binary, model4_traits_binary, model5_traits_binary, model6_traits_binary) %>% 
  map_df(tidy) %>% 
  filter(term == "MunichAnalogy" | term == "CubanMissileAnalogy" | term =="BosniaAnalogy" |
           term == "VietnamAnalogy" | term == "PhilippinesAnalogy" | term == "MunichAnalogy2")

#Meta Effect = 0.1417 and p-value = <0.0001
meta_traits <- rma.uni(yi = estimate, sei = std.error, method = "REML", data = meta_df)


##Moral Obligation
meta_df <- 
  list(model1_moral_binary, model2_moral_binary, model3_moral_binary, model4_moral_binary, model5_moral_binary, model6_moral_binary) %>% 
  map_df(tidy) %>% 
  filter(term == "MunichAnalogy" | term == "CubanMissileAnalogy" | term =="BosniaAnalogy" |
           term == "VietnamAnalogy" | term == "PhilippinesAnalogy" | term == "MunichAnalogy2")

#Meta Effect = 0.0463 and p-value = 0.0242
meta_moral <- rma.uni(yi = estimate, sei = std.error, method = "REML", data = meta_df)

######Plot Results for Analogies vs Other Types and When More Effective 
Study3$analogy_vs_gut <- NA 
Study3$analogy_vs_gut[Study3$Treatment=="Treatment Intuition"] <- 0
Study3$analogy_vs_gut[Study3$Treatment=="Treatment Munich"] <- 1

analogies_vs_gut <- lm_robust(Cred ~ analogy_vs_gut, data=Study3)

Study3$analogy_vs_expert <- NA 
Study3$analogy_vs_expert[Study3$Treatment=="Treatment Expert"] <- 0
Study3$analogy_vs_expert[Study3$Treatment=="Treatment Munich"] <- 1

analogies_vs_expert <- lm_robust(Cred ~ analogy_vs_expert, data=Study3)

Study2$pos_vs_neg_analogy <- NA
Study2$pos_vs_neg_analogy[Study2$Z1==1] <- 0
Study2$pos_vs_neg_analogy[Study2$Z1==2] <- 1

pos_vs_neg_analogy <- lm_robust(Cred_Val ~ pos_vs_neg_analogy, data=Study2)

Study2$high_vs_low_fam <- NA
Study2$high_vs_low_fam[Study2$PhilippinesAnalogy==1] <- 0
Study2$high_vs_low_fam[Study2$VietnamAnalogy==1] <- 1

high_vs_low_fam <- lm_robust(Cred_Fam ~ high_vs_low_fam, data=Study2)


estimates <- c(
  analogies_vs_gut$coefficients[-1], 
  analogies_vs_expert$coefficients[-1], 
  pos_vs_neg_analogy$coefficients[-1], 
  high_vs_low_fam$coefficients[-1]
)

conf_low <- c(
  analogies_vs_gut$conf.low[-1], 
  analogies_vs_expert$conf.low[-1], 
  pos_vs_neg_analogy$conf.low[-1], 
  high_vs_low_fam$conf.low[-1]
)

conf_high <- c(
  analogies_vs_gut$conf.high[-1], 
  analogies_vs_expert$conf.high[-1], 
  pos_vs_neg_analogy$conf.high[-1], 
  high_vs_low_fam$conf.high[-1]
)

cat("Length of estimates: ", length(estimates), "\n")
cat("Length of conf.low: ", length(conf_low), "\n")
cat("Length of conf.high: ", length(conf_high), "\n")
cat("Length of category: ", length(c(
  rep("Historical Analogies vs. Leader's Gut", length(analogies_vs_gut$coefficients[-1])),
  rep("Historical Analogies vs. Appeal to Experts", length(analogies_vs_expert$coefficients[-1])),
  rep("Positive vs. Negative Analogies", length(pos_vs_neg_analogy$coefficients[-1])),
  rep("High vs. Low Familiarity Analogies", length(high_vs_low_fam$coefficients[-1]))
)), "\n")


category <- c(
  rep("Historical Analogies vs. Leader's Gut", length(analogies_vs_gut$coefficients[-1])),
  rep("Historical Analogies vs. Appeal to Experts", length(analogies_vs_expert$coefficients[-1])),
  rep("Positive vs. Negative Analogies", length(pos_vs_neg_analogy$coefficients[-1])),
  rep("High vs. Low Familiarity Analogies", length(high_vs_low_fam$coefficients[-1]))
)

# Create the data frame
results <- data.frame(
  estimate = estimates,
  conf.low = conf_low,
  conf.high = conf_high,
  category = category
)


# Order the categories for plotting
results$category <- factor(results$category, levels = c("High vs. Low Familiarity Analogies", "Positive vs. Negative Analogies", "Historical Analogies vs. Appeal to Experts", "Historical Analogies vs. Leader's Gut"))

plot <- ggplot(results, aes(x = estimate, y = category, color = category)) +
  geom_point(size = 4) +
  geom_linerange(aes(xmin = conf.low, xmax = conf.high), size = 1) +
  theme_minimal(base_size = 15) +
  ylab("") +
  xlab("\n Effect of Historical Analogies on \n Perceptions the President Chose the \n Best Foreign Policy Strategy") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red", size = 1) +
  #geom_hline(yintercept = 2.5, linetype = "dotted", color = "black", size = 1) +
  theme(
    axis.text.x = element_text(size = 15, color = "black"),
    axis.text.y = element_text(size = 15, color = "black"),
    axis.title.y = element_text(size = 18, face = "bold"),
    axis.title.x = element_text(size = 18, face = "bold"),
    plot.caption = element_text(size = 15)
  ) +
  labs(
    title = "",
  ) +
  scale_color_manual(values = c("Historical Analogies vs. Leader's Gut" = "black", "Historical Analogies vs. Appeal to Experts" = "black", "Positive vs. Negative Analogies" = "black", "High vs. Low Familiarity Analogies" = "black")) +
  theme(legend.position = "none")

# Display the plot
print(plot)


### Heterogeneous Effects -- Study 1

Study1$co_partisans <- NA
# Replace with 0 when voter and president are from different parties
Study1$co_partisans[(Study1$democrat == 1 & Study1$republicanpresident == 1) |
                  (Study1$republican == 1 & Study1$democraticpresident == 1)] <- 0

# Replace with 1 when voter and president are from same party
Study1$co_partisans[(Study1$democrat == 1 & Study1$democraticpresident == 1) |
                  (Study1$republican == 1 & Study1$republicanpresident == 1)] <- 1

Study1$over_70 <- 0
Study1$over_70[Study1$age >= 70] <- 1

Study1$analogy_treatment <- NA 
Study1$analogy_treatment[Study1$CubanMissileAnalogy==0] <- 0
Study1$analogy_treatment[Study1$MunichAnalogy==0] <- 0
Study1$analogy_treatment[Study1$CubanMissileAnalogy==1] <- 1
Study1$analogy_treatment[Study1$MunichAnalogy==1] <- 1

# formula 
formula_pred <- bf(
  Cred ~ lambda*analogy_treatment + controls,
  
  lambda ~ Disposition_Pooled*PID*FP_Knowledge*education*Male*age*over_70, 
  
  controls ~  RepublicanPresident +Inc+ Inc_RF+ White,
  
  nl = TRUE
)

pred_prior <- c(
  prior(normal(0, 1), nlpar = "lambda"),
  prior(normal(0, .5), nlpar = "controls")
)

tw.treat.het <- brm(formula_pred,
  data = Study1,
  prior = pred_prior,
  family = gaussian(),
  cores = 4,
  # control = list(adapt_delta = .99,
  #                max_treedepth = 20),
  backend = "cmdstanr",
  refresh = 500
)
summary(tw.treat.het)

grid.treat.het <- Study1 %>%
  select(analogy_treatment, Disposition_Pooled, PID, FP_Knowledge, 
         education, Male, age, over_70) %>%
  distinct() %>%
  mutate()

# predicted outcomes
pred.treat.het <- tw.treat.het %>%
  add_epred_draws(newdata = Study1) %>%  # This adds posterior draws
  group_by(ResponseId) %>%
  select(ResponseId, .epred, analogy_treatment) %>%
  summarise(across(everything(), list(pred_median = median)), .groups = "drop")


ggplot(pred.treat.het, aes(x = .epred_pred_median, y = ResponseId, 
                           color = factor(analogy_treatment_pred_median))) +
  geom_point() +
  labs(x = "Predicted Credibility",
       y = "",
       color = "Analogy")

# Look at slopes
lambda.est <- posterior_epred(tw.treat.het, nlpar = "lambda")
lambda.quant <- apply(lambda.est, 2, 
                      function(x) quantile(x, probs = c(.1, .5, .9)))
rownames(lambda.quant) <- c("treat.10", "treat.med", "treat.90")

tw.est <- bind_cols(Study1, t(lambda.quant)) %>%
  filter(analogy_treatment == 1)
glimpse(tw.est)

ggplot(tw.est, aes(x = treat.med, y = ResponseId)) +
  geom_point() +
  labs(x = "Estimated Treatment Effect",
       y = "",
       color = "Analogy")

ggplot(tw.est, aes(x = treat.med)) +
  geom_density()

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    Disposition_Pooled
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      variable == "Disposition_Pooled" ~ "Hawkishness",
#      variable == "PID" ~ "Partisanship",
#      variable == "education" ~ "Education",
#      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
#      variable == "Male" ~ "Male",
#      variable == "age" ~ "Age",
#      variable == "over_70" ~ "Over 70"
    )
  )

hawk_study1<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(
    breaks = c("1", "2", "3", "4", "5"),
    labels = c("Very Dovish", "Dovish", "Moderate", "Hawkish", "Very Hawkish")
  )

print(hawk_study1)
ggsave("figures/hawkish_heterogeneity_study1.png", plot = hawk_study1, width = 10, height = 7, dpi = 300, device = "png")

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    PID
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      #      variable == "Disposition_Pooled" ~ "Hawkishness",
            variable == "PID" ~ "Partisanship",
      #      variable == "education" ~ "Education",
      #      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      #      variable == "Male" ~ "Male",
      #      variable == "age" ~ "Age",
      #      variable == "over_70" ~ "Over 70"
    )
  )

pid_study1<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(
    breaks = c("1", "2", "3", "4", "5", "6", "7"),
    labels = c("Strong Democrat", "Democrat", "Lean Democrat", "Independent/Other", "Lean Republican", "Republican", "Strong Republican")
  )

print(pid_study1)
ggsave("figures/pid_heterogeneity_study1.png", plot = pid_study1, width = 10, height = 7, dpi = 300, device = "png")

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    education
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      #      variable == "Disposition_Pooled" ~ "Hawkishness",
      #      variable == "PID" ~ "Partisanship",
            variable == "education" ~ "Education",
      #      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      #      variable == "Male" ~ "Male",
      #      variable == "age" ~ "Age",
      #      variable == "over_70" ~ "Over 70"
    )
  )

edu_study1<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(
    breaks = c("1", "2", "3", "4", "5", "6", "7", "8"),
    labels = c("Some High School", "High School Degree", "Vocational Degree", "Some College", "Associate's Degree", "Bachelor's Degree", "Professional Degree", "Doctoral Degree")
  )

print(edu_study1)
ggsave("figures/edu_heterogeneity_study1.png", plot = edu_study1, width = 10, height = 7, dpi = 300, device = "png")

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    FP_Knowledge
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      #      variable == "Disposition_Pooled" ~ "Hawkishness",
      #      variable == "PID" ~ "Partisanship",
      #      variable == "education" ~ "Education",
            variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      #      variable == "Male" ~ "Male",
      #      variable == "age" ~ "Age",
      #      variable == "over_70" ~ "Over 70"
    )
  )

fp_study1<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(
    breaks = c("0", "0.2", "0.4", "0.6", "0.8", "1"),
    labels = c("Very Low Knowledge", "Low Knowledge", "Some Knowledge", "Moderate Knowledge", "High Knowledge", "Very High Knowledge")
  )

print(fp_study1)
ggsave("figures/fp_heterogeneity_study1.png", plot = fp_study1, width = 10, height = 7, dpi = 300, device = "png")

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    Male
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      #      variable == "Disposition_Pooled" ~ "Hawkishness",
      #      variable == "PID" ~ "Partisanship",
      #      variable == "education" ~ "Education",
      #      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
            variable == "Male" ~ "Male",
      #      variable == "age" ~ "Age",
      #      variable == "over_70" ~ "Over 70"
    )
  )

male_study1<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(
    breaks = c("0", "1"),
    labels = c("Woman", "Man")
  )

print(male_study1)
ggsave("figures/male_heterogeneity_study1.png", plot = male_study1, width = 10, height = 7, dpi = 300, device = "png")

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    age
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      #      variable == "Disposition_Pooled" ~ "Hawkishness",
      #      variable == "PID" ~ "Partisanship",
      #      variable == "education" ~ "Education",
      #      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      #      variable == "Male" ~ "Male",
            variable == "age" ~ "Age",
      #      variable == "over_70" ~ "Over 70"
    )
  )

age_study1<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(
    breaks = c("18", "28", "38", "48", "58", "68", "78", "88"),
    labels = c("18", "28", "38", "48", "58", "68", "78", "88")
  )

print(age_study1)
ggsave("figures/age_heterogeneity_study1.png", plot = age_study1, width = 10, height = 7, dpi = 300, device = "png")

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    over_70
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      #      variable == "Disposition_Pooled" ~ "Hawkishness",
      #      variable == "PID" ~ "Partisanship",
      #      variable == "education" ~ "Education",
      #      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      #      variable == "Male" ~ "Male",
      #      variable == "age" ~ "Age",
            variable == "over_70" ~ "Aged 70+"
    )
  )

over70_study1<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(
    breaks = c("0", "1"),
    labels = c("Younger than 70", "70 or Older")
  )

print(over70_study1)
ggsave("figures/over70_heterogeneity_study1.png", plot = over70_study1, width = 10, height = 7, dpi = 300, device = "png")

### Heterogeneous Effects -- Study 3

Study3 <- Study3 %>%
  filter(!is.na(Cred))

Study3$row_id <- 1:nrow(Study3)

Study3$over_70 <- 0
Study3$over_70[Study3$Age >= 70] <- 1

Study3$analogy_treatment <- NA 
Study3$analogy_treatment[Study3$Treatment=="Control Taiwan"] <- 0
Study3$analogy_treatment[Study3$Treatment=="Treatment Expert"] <- 0
Study3$analogy_treatment[Study3$Treatment=="Treatment Intuition"] <- 0
Study3$analogy_treatment[Study3$Treatment=="Treatment Munich"] <- 1

# formula 
formula_pred <- bf(
  Cred ~ lambda*analogy_treatment + controls,
  
  lambda ~ Disposition_Pooled*PID*FP_Knowledge*College*Male*Age*over_70, 
  
  controls ~  RepublicanPresident +Inc+ White,
  
  nl = TRUE
)

pred_prior <- c(
  prior(normal(0, 1), nlpar = "lambda"),
  prior(normal(0, .5), nlpar = "controls")
)

tw.treat.het <- brm(formula_pred,
                    data = Study3,
                    prior = pred_prior,
                    family = gaussian(),
                    cores = 4,
                    # control = list(adapt_delta = .99,
                    #                max_treedepth = 20),
                    backend = "cmdstanr",
                    refresh = 500
)
summary(tw.treat.het)

grid.treat.het <- Study3 %>%
  select(analogy_treatment, Disposition_Pooled, PID, FP_Knowledge, 
         College, Male, Age, over_70) %>%
  distinct() %>%
  mutate()

# predicted outcomes
pred.treat.het <- tw.treat.het %>%
  add_epred_draws(newdata = Study3) %>%  # This adds posterior draws
  group_by(row_id) %>%
  select(row_id, .epred, analogy_treatment) %>%
  summarise(across(everything(), list(pred_median = median)), .groups = "drop")

ggplot(pred.treat.het, aes(x = .epred_pred_median, y = row_id, 
                           color = factor(analogy_treatment_pred_median))) +
  geom_point() +
  labs(x = "Predicted Credibility",
       y = "",
       color = "Analogy")

# Look at slopes
lambda.est <- posterior_epred(tw.treat.het, nlpar = "lambda")
lambda.quant <- apply(lambda.est, 2, 
                      function(x) quantile(x, probs = c(.1, .5, .9)))
rownames(lambda.quant) <- c("treat.10", "treat.med", "treat.90")

tw.est <- bind_cols(Study3, t(lambda.quant)) %>%
  filter(analogy_treatment == 1)
glimpse(tw.est)

ggplot(tw.est, aes(x = treat.med, y = row_id)) +
  geom_point() +
  labs(x = "Estimated Treatment Effect",
       y = "",
       color = "Analogy")

ggplot(tw.est, aes(x = treat.med)) +
  geom_density()

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    Disposition_Pooled
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      variable == "Disposition_Pooled" ~ "Hawkishness",
      #      variable == "PID" ~ "Partisanship",
      #      variable == "education" ~ "Education",
      #      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      #      variable == "Male" ~ "Male",
      #      variable == "Age" ~ "Age",
      #      variable == "over_70" ~ "Over 70"
    )
  )

hawk_study3<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(
    breaks = c("1", "2", "3", "4", "5"),
    labels = c("Very Dovish", "Dovish", "Moderate", "Hawkish", "Very Hawkish")
  )

print(hawk_study3)
ggsave("figures/hawkish_heterogeneity_study3.png", plot = hawk_study3, width = 10, height = 7, dpi = 300, device = "png")

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    PID
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      #      variable == "Disposition_Pooled" ~ "Hawkishness",
      variable == "PID" ~ "Partisanship",
      #      variable == "education" ~ "Education",
      #      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      #      variable == "Male" ~ "Male",
      #      variable == "Age" ~ "Age",
      #      variable == "over_70" ~ "Over 70"
    )
  )

pid_study3<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(
    breaks = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10"),
    labels = c("Strong Democrat", "Democrat", "Independent\nDemocrat", "Independent", "Independent\nRepublican", "Other,\nLean Democrat", "Other", "Other,\nLean Republican", "Republican", "Strong Republican")
  )

print(pid_study3)
ggsave("figures/pid_heterogeneity_study3.png", plot = pid_study3, width = 10, height = 7, dpi = 300, device = "png")

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    College
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      #      variable == "Disposition_Pooled" ~ "Hawkishness",
      #      variable == "PID" ~ "Partisanship",
      variable == "College" ~ "College Educated",
      #      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      #      variable == "Male" ~ "Male",
      #      variable == "Age" ~ "Age",
      #      variable == "over_70" ~ "Over 70"
    )
  )

edu_study3<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(
    breaks = c("0", "1"),
    labels = c("No College", "College Educated")
  )

print(edu_study3)
ggsave("figures/edu_heterogeneity_study3.png", plot = edu_study3, width = 10, height = 7, dpi = 300, device = "png")

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    FP_Knowledge
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      #      variable == "Disposition_Pooled" ~ "Hawkishness",
      #      variable == "PID" ~ "Partisanship",
      #      variable == "education" ~ "Education",
      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      #      variable == "Male" ~ "Male",
      #      variable == "Age" ~ "Age",
      #      variable == "over_70" ~ "Over 70"
    )
  )

fp_study3<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(
    breaks = c("0", "0.2", "0.4", "0.6", "0.8", "1"),
    labels = c("Very Low Knowledge", "Low Knowledge", "Some Knowledge", "Moderate Knowledge", "High Knowledge", "Very High Knowledge")
  )

print(fp_study3)
ggsave("figures/fp_heterogeneity_study3.png", plot = fp_study3, width = 10, height = 7, dpi = 300, device = "png")

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    Male
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      #      variable == "Disposition_Pooled" ~ "Hawkishness",
      #      variable == "PID" ~ "Partisanship",
      #      variable == "education" ~ "Education",
      #      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      variable == "Male" ~ "Male",
      #      variable == "Age" ~ "Age",
      #      variable == "over_70" ~ "Over 70"
    )
  )

male_study3<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(
    breaks = c("0", "1"),
    labels = c("Woman", "Man")
  )

print(male_study3)
ggsave("figures/male_heterogeneity_study3.png", plot = male_study3, width = 10, height = 7, dpi = 300, device = "png")

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    Age
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      #      variable == "Disposition_Pooled" ~ "Hawkishness",
      #      variable == "PID" ~ "Partisanship",
      #      variable == "education" ~ "Education",
      #      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      #      variable == "Male" ~ "Male",
      variable == "Age" ~ "Age",
      #      variable == "over_70" ~ "Over 70"
    )
  )

age_study3<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(
    breaks = c("18", "28", "38", "48", "58", "68", "78", "88"),
    labels = c("18", "28", "38", "48", "58", "68", "78", "88")
  )

print(age_study3)
ggsave("figures/age_heterogeneity_study3.png", plot = age_study3, width = 10, height = 7, dpi = 300, device = "png")

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    over_70
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      #      variable == "Disposition_Pooled" ~ "Hawkishness",
      #      variable == "PID" ~ "Partisanship",
      #      variable == "education" ~ "Education",
      #      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      #      variable == "Male" ~ "Male",
      #      variable == "Age" ~ "Age",
      variable == "over_70" ~ "Aged 70+"
    )
  )

over70_study3<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(
    breaks = c("0", "1"),
    labels = c("Younger than 70", "70 or Older")
  )

print(over70_study3)
ggsave("figures/over70_heterogeneity_study3.png", plot = over70_study3, width = 10, height = 7, dpi = 300, device = "png")


### Heterogeneous Effects -- Study 2 -- Myanmar

Study2$over_40 <- 0
Study2$over_40[Study2$Age >= 40] <- 1

Study2$analogy_treatment <- NA 
Study2$analogy_treatment[Study2$Z1=="0"] <- 0
Study2$analogy_treatment[Study2$Z1=="1"] <- 1
Study2$analogy_treatment[Study2$Z1=="2"] <- 1
Study2$analogy_treatment[Study2$Z1=="3"] <- 0
Study2$analogy_treatment[Study2$Z1=="4"] <- 1
Study2$analogy_treatment[Study2$Z1=="5"] <- 1

Study2$myanmar <- NA 
Study2$myanmar[Study2$Z1=="0"] <- 1
Study2$myanmar[Study2$Z1=="1"] <- 1
Study2$myanmar[Study2$Z1=="2"] <- 1
Study2$myanmar[Study2$Z1=="3"] <- 0
Study2$myanmar[Study2$Z1=="4"] <- 0
Study2$myanmar[Study2$Z1=="5"] <- 0

Study2Myanmar <- Study2 %>%
  filter(myanmar==1)

Study2Myanmar$Cred_Val <- as.numeric(Study2Myanmar$Cred_Val)

# formula 
formula_pred <- bf(
  Cred_Val ~ lambda*analogy_treatment + controls,
  
  lambda ~ Negativity_Bias_Pooled*Disposition_Pooled*PID*FP_Knowledge*education*Male*Age*over_40, 
  
  controls ~  RepublicanPresident +Income+ White,
  
  nl = TRUE
)

pred_prior <- c(
  prior(normal(0, 1), nlpar = "lambda"),
  prior(normal(0, .5), nlpar = "controls")
)

tw.treat.het <- brm(formula_pred,
                    data = Study2Myanmar,
                    prior = pred_prior,
                    family = gaussian(),
                    cores = 4,
                    # control = list(adapt_delta = .99,
                    #                max_treedepth = 20),
                    backend = "cmdstanr",
                    refresh = 500
)
summary(tw.treat.het)

grid.treat.het <- Study2Myanmar %>%
  select(analogy_treatment, Negativity_Bias_Pooled, Disposition_Pooled,
         PID, FP_Knowledge, education, Male, Age, over_40) %>%
  distinct() %>%
  mutate()

# predicted outcomes
pred.treat.het <- tw.treat.het %>%
  add_epred_draws(newdata = Study2Myanmar) %>%  # This adds posterior draws
  group_by(ResponseId) %>%
  select(ResponseId, .epred, analogy_treatment) %>%
  summarise(across(everything(), list(pred_median = median)), .groups = "drop")

ggplot(pred.treat.het, aes(x = .epred_pred_median, y = ResponseId, 
                           color = factor(analogy_treatment_pred_median))) +
  geom_point() +
  labs(x = "Predicted Credibility",
       y = "",
       color = "Analogy")

# Look at slopes
lambda.est <- posterior_epred(tw.treat.het, nlpar = "lambda")
lambda.quant <- apply(lambda.est, 2, 
                      function(x) quantile(x, probs = c(.1, .5, .9)))
rownames(lambda.quant) <- c("treat.10", "treat.med", "treat.90")

tw.est <- bind_cols(Study2Myanmar, t(lambda.quant)) %>%
  filter(analogy_treatment == 1)
glimpse(tw.est)

ggplot(tw.est, aes(x = treat.med, y = ResponseId)) +
  geom_point() +
  labs(x = "Estimated Treatment Effect",
       y = "",
       color = "Analogy")

ggplot(tw.est, aes(x = treat.med)) +
  geom_density()

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    Disposition_Pooled
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      variable == "Disposition_Pooled" ~ "Hawkishness",
      #      variable == "PID" ~ "Partisanship",
      #      variable == "education" ~ "Education",
      #      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      #      variable == "Male" ~ "Male",
      #      variable == "Age" ~ "Age",
      #      variable == "over_40" ~ "Over 40"
    )
  )

hawk_study2m<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(
    breaks = c("1", "2", "3", "4", "5"),
    labels = c("Very Dovish", "Dovish", "Moderate", "Hawkish", "Very Hawkish")
  )

print(hawk_study2m)
ggsave("figures/hawkish_heterogeneity_study2m.png", plot = hawk_study2m, width = 10, height = 7, dpi = 300, device = "png")

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    PID
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      #      variable == "Disposition_Pooled" ~ "Hawkishness",
      variable == "PID" ~ "Partisanship",
      #      variable == "education" ~ "Education",
      #      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      #      variable == "Male" ~ "Male",
      #      variable == "Age" ~ "Age",
      #      variable == "over_70" ~ "Over 70"
    )
  )

pid_study2m<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(
    breaks = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10"),
    labels = c("Strong Democrat", "Democrat", "Independent\nDemocrat", "Independent", "Independent\nRepublican", "Other,\nLean Democrat", "Other", "Other,\nLean Republican", "Republican", "Strong Republican")
  )

print(pid_study2m)
ggsave("figures/pid_heterogeneity_study2m.png", plot = pid_study2m, width = 10, height = 7, dpi = 300, device = "png")

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    education
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      #      variable == "Disposition_Pooled" ~ "Hawkishness",
      #      variable == "PID" ~ "Partisanship",
      variable == "education" ~ "Education",
      #      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      #      variable == "Male" ~ "Male",
      #      variable == "Age" ~ "Age",
      #      variable == "over_70" ~ "Over 70"
    )
  )

edu_study2m<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(breaks = c("1", "2", "3", "4", "5", "6", "7", "8"),
  labels = c("Some High School", "High School Degree", "Vocational Degree", "Some College", "Associate's Degree", "Bachelor's Degree", "Professional Degree", "Doctoral Degree")
      )

print(edu_study2m)
ggsave("figures/edu_heterogeneity_study2m.png", plot = edu_study2m, width = 10, height = 7, dpi = 300, device = "png")

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    FP_Knowledge
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      #      variable == "Disposition_Pooled" ~ "Hawkishness",
      #      variable == "PID" ~ "Partisanship",
      #      variable == "education" ~ "Education",
      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      #      variable == "Male" ~ "Male",
      #      variable == "Age" ~ "Age",
      #      variable == "over_70" ~ "Over 70"
    )
  )

fp_study2m<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(
    breaks = c("0", "0.2", "0.4", "0.6", "0.8", "1"),
    labels = c("Very Low Knowledge", "Low Knowledge", "Some Knowledge", "Moderate Knowledge", "High Knowledge", "Very High Knowledge")
  )

print(fp_study2m)
ggsave("figures/fp_heterogeneity_study2m.png", plot = fp_study2m, width = 10, height = 7, dpi = 300, device = "png")

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    Male
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      #      variable == "Disposition_Pooled" ~ "Hawkishness",
      #      variable == "PID" ~ "Partisanship",
      #      variable == "education" ~ "Education",
      #      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      variable == "Male" ~ "Male",
      #      variable == "Age" ~ "Age",
      #      variable == "over_70" ~ "Over 70"
    )
  )

male_study2m<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(
    breaks = c("0", "1"),
    labels = c("Woman", "Man")
  )

print(male_study2m)
ggsave("figures/male_heterogeneity_study2m.png", plot = male_study2m, width = 10, height = 7, dpi = 300, device = "png")

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    Age
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      #      variable == "Disposition_Pooled" ~ "Hawkishness",
      #      variable == "PID" ~ "Partisanship",
      #      variable == "education" ~ "Education",
      #      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      #      variable == "Male" ~ "Male",
      variable == "Age" ~ "Age",
      #      variable == "over_70" ~ "Over 70"
    )
  )

age_study2m<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(
    breaks = c("18", "28", "38", "48", "58", "68", "78", "88"),
    labels = c("18", "28", "38", "48", "58", "68", "78", "88")
  )

print(age_study2m)
ggsave("figures/age_heterogeneity_study2m.png", plot = age_study2m, width = 10, height = 7, dpi = 300, device = "png")

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    over_40
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      #      variable == "Disposition_Pooled" ~ "Hawkishness",
      #      variable == "PID" ~ "Partisanship",
      #      variable == "education" ~ "Education",
      #      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      #      variable == "Male" ~ "Male",
      #      variable == "Age" ~ "Age",
      variable == "over_40" ~ "Aged 40+"
    )
  )

over40_study2m<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(
    breaks = c("0", "1"),
    labels = c("Younger than 40", "40 or Older")
  )

print(over40_study2m)
ggsave("figures/over40_heterogeneity_study2m.png", plot = over40_study2m, width = 10, height = 7, dpi = 300, device = "png")

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    Negativity_Bias_Pooled
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      #      variable == "Disposition_Pooled" ~ "Hawkishness",
      #      variable == "PID" ~ "Partisanship",
      #      variable == "education" ~ "Education",
      #      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      #      variable == "Male" ~ "Male",
      #      variable == "Age" ~ "Age",
      #variable == "over_40" ~ "Aged 40+",
      variable == "Negativity_Bias_Pooled" ~ "Negativity Bias"
    )
  )

negativity_study2m<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  )

print(negativity_study2m)
ggsave("figures/negativity_heterogeneity_study2m.png", plot = negativity_study2m, width = 10, height = 7, dpi = 300, device = "png")

### Heterogeneous Effects -- Study 2 -- Iran

Study2$iran <- NA 
Study2$iran[Study2$Z1=="0"] <- 0
Study2$iran[Study2$Z1=="1"] <- 0
Study2$iran[Study2$Z1=="2"] <- 0
Study2$iran[Study2$Z1=="3"] <- 1
Study2$iran[Study2$Z1=="4"] <- 1
Study2$iran[Study2$Z1=="5"] <- 1

Study2Iran <- Study2 %>%
  filter(iran==1)

Study2Iran$Cred_Fam <- as.numeric(Study2Iran$Cred_Fam)

Study2Iran <- Study2Iran %>%
  filter(!is.na(Cred_Fam))

Study2Iran$over_70 <- 0
Study2Iran$over_70[Study2Iran$Age >= 70] <- 1

# formula 
formula_pred <- bf(
  Cred_Fam ~ lambda*analogy_treatment + controls,
  
  lambda ~ Negativity_Bias_Pooled*Disposition_Pooled*PID*FP_Knowledge*education*Male*Age*over_70, 
  
  controls ~  RepublicanPresident +Income+ White,
  
  nl = TRUE
)

pred_prior <- c(
  prior(normal(0, 1), nlpar = "lambda"),
  prior(normal(0, .5), nlpar = "controls")
)

tw.treat.het <- brm(formula_pred,
                    data = Study2Iran,
                    prior = pred_prior,
                    family = gaussian(),
                    cores = 4,
                    # control = list(adapt_delta = .99,
                    #                max_treedepth = 20),
                    backend = "cmdstanr",
                    refresh = 500
)
summary(tw.treat.het)

grid.treat.het <- Study2Iran %>%
  select(analogy_treatment, Negativity_Bias_Pooled, Disposition_Pooled,
         PID, FP_Knowledge, education, Male, Age, over_70) %>%
  distinct() %>%
  mutate()

# predicted outcomes
pred.treat.het <- tw.treat.het %>%
  add_epred_draws(newdata = Study2Iran) %>%  # This adds posterior draws
  group_by(ResponseId) %>%
  select(ResponseId, .epred, analogy_treatment) %>%
  summarise(across(everything(), list(pred_median = median)), .groups = "drop")

ggplot(pred.treat.het, aes(x = .epred_pred_median, y = ResponseId, 
                           color = factor(analogy_treatment_pred_median))) +
  geom_point() +
  labs(x = "Predicted Credibility",
       y = "",
       color = "Analogy")

# Look at slopes
lambda.est <- posterior_epred(tw.treat.het, nlpar = "lambda")
lambda.quant <- apply(lambda.est, 2, 
                      function(x) quantile(x, probs = c(.1, .5, .9)))
rownames(lambda.quant) <- c("treat.10", "treat.med", "treat.90")

tw.est <- bind_cols(Study2Iran, t(lambda.quant)) %>%
  filter(analogy_treatment == 1)
glimpse(tw.est)

ggplot(tw.est, aes(x = treat.med, y = ResponseId)) +
  geom_point() +
  labs(x = "Estimated Treatment Effect",
       y = "",
       color = "Analogy")

ggplot(tw.est, aes(x = treat.med)) +
  geom_density()

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    Disposition_Pooled
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      variable == "Disposition_Pooled" ~ "Hawkishness",
      #      variable == "PID" ~ "Partisanship",
      #      variable == "education" ~ "Education",
      #      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      #      variable == "Male" ~ "Male",
      #      variable == "Age" ~ "Age",
      #      variable == "over_40" ~ "Over 40"
    )
  )

hawk_study2i<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(
    breaks = c("1", "2", "3", "4", "5"),
    labels = c("Very Dovish", "Dovish", "Moderate", "Hawkish", "Very Hawkish")
  )

print(hawk_study2i)
ggsave("figures/hawkish_heterogeneity_study2i.png", plot = hawk_study2i, width = 10, height = 7, dpi = 300, device = "png")

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    PID
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      #      variable == "Disposition_Pooled" ~ "Hawkishness",
      variable == "PID" ~ "Partisanship",
      #      variable == "education" ~ "Education",
      #      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      #      variable == "Male" ~ "Male",
      #      variable == "Age" ~ "Age",
      #      variable == "over_70" ~ "Over 70"
    )
  )

pid_study2i<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(
    breaks = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10"),
    labels = c("Strong Democrat", "Democrat", "Independent\nDemocrat", "Independent", "Independent\nRepublican", "Other,\nLean Democrat", "Other", "Other,\nLean Republican", "Republican", "Strong Republican")
  )

print(pid_study2i)
ggsave("figures/pid_heterogeneity_study2i.png", plot = pid_study2i, width = 10, height = 7, dpi = 300, device = "png")

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    education
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      #      variable == "Disposition_Pooled" ~ "Hawkishness",
      #      variable == "PID" ~ "Partisanship",
      variable == "education" ~ "Education",
      #      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      #      variable == "Male" ~ "Male",
      #      variable == "Age" ~ "Age",
      #      variable == "over_70" ~ "Over 70"
    )
  )

edu_study2i<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(breaks = c("1", "2", "3", "4", "5", "6", "7", "8"),
                   labels = c("Some High School", "High School Degree", "Vocational Degree", "Some College", "Associate's Degree", "Bachelor's Degree", "Professional Degree", "Doctoral Degree")
  )

print(edu_study2i)
ggsave("figures/edu_heterogeneity_study2i.png", plot = edu_study2i, width = 10, height = 7, dpi = 300, device = "png")

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    FP_Knowledge
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      #      variable == "Disposition_Pooled" ~ "Hawkishness",
      #      variable == "PID" ~ "Partisanship",
      #      variable == "education" ~ "Education",
      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      #      variable == "Male" ~ "Male",
      #      variable == "Age" ~ "Age",
      #      variable == "over_70" ~ "Over 70"
    )
  )

fp_study2i<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(
    breaks = c("0", "0.2", "0.4", "0.6", "0.8", "1"),
    labels = c("Very Low Knowledge", "Low Knowledge", "Some Knowledge", "Moderate Knowledge", "High Knowledge", "Very High Knowledge")
  )

print(fp_study2i)
ggsave("figures/fp_heterogeneity_study2i.png", plot = fp_study2i, width = 10, height = 7, dpi = 300, device = "png")

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    Male
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      #      variable == "Disposition_Pooled" ~ "Hawkishness",
      #      variable == "PID" ~ "Partisanship",
      #      variable == "education" ~ "Education",
      #      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      variable == "Male" ~ "Male",
      #      variable == "Age" ~ "Age",
      #      variable == "over_70" ~ "Over 70"
    )
  )

male_study2i<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(
    breaks = c("0", "1"),
    labels = c("Woman", "Man")
  )

print(male_study2i)
ggsave("figures/male_heterogeneity_study2i.png", plot = male_study2i, width = 10, height = 7, dpi = 300, device = "png")

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    Age
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      #      variable == "Disposition_Pooled" ~ "Hawkishness",
      #      variable == "PID" ~ "Partisanship",
      #      variable == "education" ~ "Education",
      #      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      #      variable == "Male" ~ "Male",
      variable == "Age" ~ "Age",
      #      variable == "over_70" ~ "Over 70"
    )
  )

age_study2i<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(
    breaks = c("18", "28", "38", "48", "58", "68", "78", "88"),
    labels = c("18", "28", "38", "48", "58", "68", "78", "88")
  )

print(age_study2i)
ggsave("figures/age_heterogeneity_study2i.png", plot = age_study2i, width = 10, height = 7, dpi = 300, device = "png")

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    over_70
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      #      variable == "Disposition_Pooled" ~ "Hawkishness",
      #      variable == "PID" ~ "Partisanship",
      #      variable == "education" ~ "Education",
      #      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      #      variable == "Male" ~ "Male",
      #      variable == "Age" ~ "Age",
      variable == "over_70" ~ "Aged 70+"
    )
  )

over70_study2i<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  ) +
  scale_x_discrete(
    breaks = c("0", "1"),
    labels = c("Younger than 70", "70 or Older")
  )

print(over70_study2i)
ggsave("figures/over70_heterogeneity_study2i.png", plot = over70_study2i, width = 10, height = 7, dpi = 300, device = "png")

### look at splits by variable
slopes.treat.het.long <- tw.est %>% 
  select(
    treat.med, treat.10, treat.90,
    Negativity_Bias_Pooled
  ) %>%
  pivot_longer(cols = -c(treat.10, treat.90,
                         treat.med),
               names_to = "variable") %>%
  mutate(
    variable = case_when(
      #      variable == "Disposition_Pooled" ~ "Hawkishness",
      #      variable == "PID" ~ "Partisanship",
      #      variable == "education" ~ "Education",
      #      variable == "FP_Knowledge" ~ "Foreign Policy Knowledge",
      #      variable == "Male" ~ "Male",
      #      variable == "Age" ~ "Age",
      #variable == "over_40" ~ "Aged 40+",
      variable == "Negativity_Bias_Pooled" ~ "Negativity Bias"
    )
  )

negativity_study2i<-ggplot(slopes.treat.het.long, aes(x = factor(value), y = treat.med)) +
  facet_wrap(~ variable, scales = "free_x") +
  geom_hline(yintercept = 0) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitter(width = .25)) +
  labs(
    y = "Effect of Historical Analogy",
    x = "Modifier Value",
    title = "Variation in Impact of Historical Analogy by Grouping Variable"
  )

print(negativity_study2i)
ggsave("figures/negativity_heterogeneity_study2i.png", plot = negativity_study2i, width = 10, height = 7, dpi = 300, device = "png")
